% Estimates fixed effects model of choice under risk and temporal delay

% based on simple_preference_estimation_parametrized_rev4

clc
clear
clear global
cd ''

global  ind_n theta_xy lottery_n lottery_choice_ind   ind est_risk r_scale a1 b1 p_a1 a2 b2 p_a2    sig_th_lower_limit time     x1 x2 t1 t2  first_stage   time_choice_n time_choice_ind 


%%%%%%%%%%%%%%%%%%%%%%%% Choice Data Import %%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%% Lotteries %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%% Choice Data %%%%%%%%%%%%%%
%%% Import "lottery_choice.out" as a "Numeric Matrix", comma-delimited.
%%% Range A2:BC1225

% in indifference_threshold_calculation, can switch to whether CRRA, CARA,
% or Expo-Power Utility is used.
CARA=0
ep=0

run indifference_threshold_calculation
theta_xy=omega_th_final1';
lottery_n=size(theta_xy,2);

% lottery one is the safer one
a1=lotteries(:,1)';
b1=lotteries(:,2)';
p_a1=lotteries(:,3)';
% lottery 2 is the riskier one
a2=lotteries(:,4)';
b2=lotteries(:,5)';
p_a2=lotteries(:,6)';

lottery_choice_ind=lotterychoice;

%%%%%%%%%%%%%%%%%%%%%%%%%% Time Choices %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%% Temporal Choices %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%% Choice Data %%%%%%%%%%%%%%
%%% Import "temporal_choice.out" as a "Numeric Matrix", comma-delimited.
%%% Range A2:AV1225

load temporal_task_info
time_choice_n=size(choices,1);

% choice 1 is the sooner one
% $ in first option
x1=choices(:,1)';
% years from now that the $ will be received
t1=choices(:,2)';
% choice 2 is the more distant one
% $ in second option
x2=choices(:,3)';
% years from now that the $ will be received
t2=choices(:,4)';

time_choice_ind=temporalchoice;

% upper limit on discount rate
r_scale=1;
% lower limit on risk aversion sd
sig_th_lower_limit=.001;


%%% number of individuals
ind_n=1224;
%%% creates an index for the individuals
ind_id = kron(ones(1,1), kron([1:ind_n]', ones(1,1))); 

%%% Additional estimation and optimization parameters 
iter_n=50;



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%Estimation%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


rng('shuffle')






%%%%%%%%%%%%%%%%%%%%%% STEP 1 - Risk only - Individual level %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%% for greater ease, put combined time estimation in one file also. The
%%% below indicator toggles whether only risk or also time is estimated.
%%% Above switches (MV, sig_exp_RU, no_E, no_sig) still work. If want a
%%% disction between risk and time sig, need to put design (risk=0, time=1)
%%% into "sig_vars"
time=0

%%% How many risk parameters to estimate 
size_risk=3;

est_risk=NaN(ind_n,size_risk);
est_risk_se=NaN(ind_n,size_risk);
exit_table1=NaN(ind_n,1);   
fval_min_table1=NaN(ind_n,1);   

theta_start_risk_only=[randn(size_risk,1)];




% display off
options = optimset('Algorithm','sqp','Display','off','TolFun',1e-200, 'TolX', 1e-200, 'MaxFunEvals', 100000, 'MaxIter', 200, 'ObjectiveLimit', -1e200);

%%% optimization individual by individual
%%% Risk estimation
for ind=1:ind_n

    result_table_all1=NaN(size_risk,iter_n);
    %%% do NaN to distinguish between failed optimization and 0
    fval_table1=NaN(1,iter_n);
    nan_table1=NaN(1,iter_n);
    start_table_all1=NaN(size_risk,iter_n);    
    i=1;
    tic
    while i<=iter_n
                    try
                        rng('shuffle')
                        
                        theta_start_risk_only=[randn(size_risk,1)];

                        [B_ML_probit_all_table, fval] = fminunc(@ml_risk_simple_ind_cdf,theta_start_risk_only, options);
                        catch ME 
                        i
                        continue;  
                    end;
          
                    
          result_table_all1(:,i)=B_ML_probit_all_table;
          start_table_all1(:,i)=theta_start_risk_only;  
                    

          fval_table1(1,i)=fval; 
          i=i+1;
    end;
    toc
    
    %%% Gets the position(s) of the cell(s) with the lowest fval
    [xxx,yyy]=find(fval_table1==min(fval_table1));
    %%% Gets the estimated coeffs with the lowest fval and highest
    %%% precision
    reduced_vector=result_table_all1(:,yyy); 
    [xxx2, yyy2]=find(reduced_vector==min(reduced_vector(3,:)));
    yyy=yyy(yyy2);
    reduced_vector=result_table_all1(:,yyy);
    [xxx2, yyy2]=find(reduced_vector==min(reduced_vector(2,:)));
    yyy=yyy(yyy2);
    reduced_vector=result_table_all1(:,yyy); 
    [xxx2,yyy2]=find(abs(reduced_vector)==min(abs(reduced_vector(1,:))));
    yyy=yyy(yyy2);
    reduced_vector=result_table_all1(:,yyy);
    col=yyy(min(yyy2));   
    step1_coeffs=result_table_all1(:,col);
    fval_min_table1(ind)=fval_table1(col);
    est_risk(ind,1:size_risk)=step1_coeffs';
 
    
    ind
    savefile = 'fe_risk.mat';
    save(savefile, 'result_table_all1', 'fval_min_table1', 'est_risk', 'iter_n', 'sig_th_lower_limit');

end    

%%% converts raw coefficients into parameters of interest
conv_risk=[est_risk(:,1) 0.5*normcdf(est_risk(:,2)) normcdf(est_risk(:,3))];
est_theta_i=conv_risk(:,1);
est_sig_th_i=conv_risk(:,3);
est_kappa_i=conv_risk(:,2);

%%% gives basic summary statistics
pct25_l=[prctile(est_theta_i, 25)  prctile(est_sig_th_i, 25) prctile(est_kappa_i, 25)];
pct50_l=[prctile(est_theta_i, 50)  prctile(est_sig_th_i, 50) prctile(est_kappa_i, 50)];
pct75_l=[prctile(est_theta_i, 75)  prctile(est_sig_th_i, 75) prctile(est_kappa_i, 75)];

pct=[pct25_l; pct50_l; pct75_l]




%%% How many time paras to estimate 
est_n_time=2;
size_time=2;
theta_start_time=[randn(size_time,1)];
est_time=NaN(ind_n,est_n_time);
exit_table2=NaN(ind_n,1);   
fval_min_table2=NaN(ind_n,1); 

%%% Time estimation
for ind=1:ind_n
    %%% risk coefficient estimated in the previous step
    first_stage=est_risk(ind,:);
  
    result_table_all2=NaN(size(theta_start_time,1),iter_n);
    fval_table2=NaN(1,iter_n);
    start_table_all2=NaN(size(theta_start_time,1),iter_n);  
    i=1;
    tic
    while i<=iter_n
                    try
                        rng('shuffle')
                        theta_start_time=[randn(size_time,1)];
                        [B_ML_probit_all_table, fval] = fminunc(@ml_time_simple_ind_cdf,theta_start_time, options);                                   
                        catch ME 
                        iter_n
                        continue; 
                    end;
          result_table_all2(:,i)=B_ML_probit_all_table;
          fval_table2(1,i)=fval;  
          i=i+1;
    end;
    toc
    
    fval_table2(fval_table2==0)=NaN;
    [M,I]=min(fval_table2);
    col=I;
    step2_coeffs=result_table_all2(:,col);
    fval_min_table2(ind)=M;
    est_time(ind,1:est_n_time)=step2_coeffs';
    

    ind

    savefile = 'fe_time.mat';
    save(savefile, 'result_table_all2', 'fval_min_table2', 'exit_table2', 'est_time');
    
    

    
end


%%% converts raw coefficients into parameters of interest
conv_time=[ (r_scale*normcdf(est_time(:,1))) normcdf(est_time(:,2))];

est_r_i=conv_time(:,1);
est_sig_r_i=conv_time(:,2);

%%% gives basic summary statistics
pct25_l=[prctile(est_theta_i, 25)  prctile(est_r_i, 25) prctile(est_sig_th_i, 25) prctile(est_sig_r_i, 25) prctile(est_kappa_i, 25)];
pct50_l=[prctile(est_theta_i, 50)  prctile(est_r_i, 50) prctile(est_sig_th_i, 50) prctile(est_sig_r_i, 50) prctile(est_kappa_i, 50)];
pct75_l=[prctile(est_theta_i, 75)  prctile(est_r_i, 75) prctile(est_sig_th_i, 75) prctile(est_sig_r_i, 75) prctile(est_kappa_i, 75)];

pct=[pct25_l; pct50_l; pct75_l]

%%% Sets bounds
est_theta_i(est_theta_i<-3)=-3;
est_theta_i(est_theta_i>5)=5;

%%% export into stata for analysis
exp_matrix=[ind_id est_theta_i est_sig_th_i est_kappa_i est_r_i est_sig_r_i fval_min_table1];
filename = 'fe_paras.csv';
csvwrite(filename,exp_matrix) 


%%% Generates moments of the choice data
time_choice=temporalchoice;
lottery_choice=lotterychoice;
run reversal_index


%%% higher values indicate more risky or more distant choices
time_index=sum(time_choice,2);
lottery_index=sum(lottery_choice,2);

%%% Export moments of the choice data
exp_matrix=[ind_id   time_index lottery_index time_reversal_n risk_reversal_n reversal_n  risk_first_ones time_first_ones ];
filename = 'choice_moments.csv';
csvwrite(filename,exp_matrix)  
